Monte Carlo Simulations of an Extended Feynman-Kikuchi Model 
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We present Quantum Monte Carlo simulations of a generalization of the Feynman-Kikuchi model 
which includes the possibility of vacancies and interactions between the particles undergoing ex- 
change. By measuring the winding number (superfluid density) and density structure factor, we 
determine the phase diagram, and show that it exhibits regions which possess both superfluid and 
charge ordering. 

PACS numbers: 05.30.Jp, 03.75.Hh, 



I. INTRODUCTION 

The study of continuum superfluid phase transitions 
using Quantum Monte Carlo (QMC) methods has a his- 
tory which includes path integral simulations of Helium 
using realistic interatomic potentials, which capture Tx 
in good quantitative agreement with experimenti, to re- 
cent numerical work^ focusing on experiments^!^ which 
observe 'supersolid' order^, the simultaneous presence of 
both superfluidity and long range density correlations. 

At the same time, related path integral studies of lat- 
tice models (the 'boson-Hubbard' Hamiltonian)^il have 
been undertaken. These too have been partially mo- 
tivated by the issue of supersolid order, but have also 
been driven by the possibility of studying the uni- 
versal conductivity in granular superconductors^', and 
superfluid-Mott insulator transitions of relevance to op- 
tically trapped atomsl 

Many of these simulations emphasize Feynman's pic- 
ture of the connection between the superfluid transition 
and the increasing entanglement (and ultimate develop- 
ment of macroscopic 'winding' across the whole sample) 
of quantum paths as the temperature is lowered. Indeed, 
the superfluid density ps is proportional to the mean 
square winding of paths around the lattice^. 

However, even before the advent of these large scale 
QMC simulations which allow the study of the superfluid 
transition exactly, Feynman^ and Kikuchii^i^^ suggested, 
and studied analytically, an approximate 'classical' model 
whose configurations are permutation loops of sites on a 
d = 3 lattice. The partition function they suggested is, 

N 



xexp[-I^5:(r.-7'r,0^]. (1) 
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Here is the effective mass of He atoms, and the func- 
tion p(ri, r2, ■ • • , tn) is assumed to be nonvanishing only 
when the coordinates are located on the sites of a reg- 
ular lattice (a cubic lattice in the original treatments). V 
refers to a permutation of the coordinates. A transition 
to a 'superfluid' phase where a macroscopic number of 
sites participate in a single large loop, as the tempera- 
ture T is lowered, was discovered and investigated. 



Various approximations were employed to determine 
the properties of this model, predominantly diagram- 
matic (series) expansions in the exchange loops. To facli- 
tate these analytic treatments, in many of the early stud- 
ies the allowed permutations were restricted to "near- 
neighbor" exchange in which the maximum "distance 
traveled" by each particle i, is only one lattice constant 
d. That is, jr^ — T^r^j = d for all particles i. One of the 
early issues concerned whether the superfluid transition 
was third order, as originally found by Feynman^., or sec- 
ond orde r^^d^d'^d^ and whether the order was affected 
by the restriction to near-neighbor exchange. Another 
issue was the behavior of the speciflc heat, both how to 
eliminate various artificial structures (and even negative 
values) near the phase transition, and also how to recover 
the experimentally observed behavior in superfluids at 
low temperatures. Here the removal of the restriction to 
local permutations was found to be crucialii. 



The Feynman-Kikuchi (FK) model is closely related to 
the duality-transformed XY modeU^iii, where the par- 
tition function can also be expressed in terms of sums 
of closed paths on a lattice. The allowed conflgurations 
are somewhat different, since in the XY case path over- 
lap is allowed whereas in the FK model, each appears 
only once in Vvi. However, the energy for paths grows 
quadratically with the overlap, so that in practice large 
overlaps are unlikely, enhancing the similarities between 
the partition functions. This connection is perhaps not 
so surprising since both models offer ways to understand 
the superfluid phase transition. 



In this paper we will study the FK model in d = 2 
using Monte Carlo simulationsil. Motivated by recent 
work on supersolids, we will then suggest a generalization 
which contains 'vacancies' and interactions between the 
occupied sites. We will determine the nature of the su- 
perfluid phase transition, and how it depends on particle 
density, and also study the possibility of charge ordered 
states arising from the interactions. The results allow 
up to construct the phase diagram of our generalized FK 
model. 
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II. MODEL AND COMPUTATIONAL 
METHODS 

We begin by briefly reviewing the motivation for the 
FK model which will expose the connection with exact 
path integral expressions for the partition function. Con- 
sider the quantum Hamiltonian for a system of N inter- 
acting bosons, 

N ^2 

^ = E|^ + ^(fl'f2,---fA^) . (2) 

i—1 

Here and f ^ are the momentum and position operators. 
The partition function is given by, 

Z = Tre-^/^ = Tr[e-^^/^... e-^^/^] 

« Tr [e-^^/^ e-'^/^ • • • e'^^/^ e"^^/^ ] (3) 

where, following the usual the path integral 
approacbi^ii^, a small parameter e has introduced, 
the exponential of the full Hamiltonian has been broken 
into M pieces with Me = 1, and then approximated 
by the product of the exponentials of the kinetic 
and potential energies individually. This 'Trotter' 
approximation32i2ii2^ becomes exact in the limit e ^ 
(A/ oo). We have set Boltzmann's constant fee = 1 
for simplicity. 

The trace is evaluated by summing over a complete set 
of position eigenstates, and also inserting additional com- 
plete sets of position eigenstates throughout the string 
of incremental imaginary time evolution operators. The 
potential energy exponentials act on the eigenstates to 
give numbers, and the remaining matrix elements of the 
kinetic energy operators are readily computed, yielding, 

» N M 

-p i=\ m=l 

M 

m— 1 

+ ^EE[^-77^]- (4) 

m— 1 i—1 

Here the superscript m is an 'imaginary time' index 
which labels the point of insertion of the different com- 
plete sets of states. The final set of positions {rf^} is 
constrained to be a permutation V of the original posi- 
tions {v}}, as a consequence of the trace in the definition 
of the quantum partition function. The sum over per- 
mutations V incorporates the indistiguishability of the 
bosonic particles. This completes the representation of 
the partition function as an integral over classical paths 
in space and imaginary time. 

Examination of this exact expression for the partition 
function, now readily motivates the origin of the FK 
model. The function /5(ri, r2, • • • , tat), and its restric- 
tion to a cubic lattice, can be thought of as arising from 



the potential energy terms which act to tend to localize 
the particle positions in a regular array. The single expo- 
nential in the particle positions in the FK model can be 
regarded as a truncation of the complete set of M expo- 
nentials for all imaginary times. Notice that the combi- 
nation of the leading factor e/T and the two such factors 
in the denominator of the 'kinetic energy' lead to the ap- 
pearance of the temperature in the numerator. From the 
viewpoint of the original path integral, this refiects the 
fact that as the temperature is lowered, the paths have 
more (imaginary) time in which to propagate, and it be- 
comes increasingly easy for them to permute. Physically, 
this then leads to a superfluid phase transition as T is 
lowered. 

In this paper, we will work with the FK model on a 
two dimensional square lattice. We will denote by (3 the 
inverse of the prefactor of the sum of the distances trav- 
eled by the individual particles. That is, our partition 
function will be, 

^ = E -p[ - f ] - E -p [ - ^ E(r^ - ^-0^] (5) 

V ' V ' i 

To be explicit, to a labeling of the sites of the two dimen- 
sional square lattice, illustrated in Fig. 1 (left) we asso- 
ciate a permutation V. An example is shown in Fig. 1 
(right) . This configuration contains three nontrivial per- 
mutation loops: The particles at sites 8 and 9 are in a 
'two particle' loop, as are the particles at sites 5 and 
18. Note that in the former each particle moves a single 
lattice site, and contributes '1' to the energy, while the 
square of the distances traveled by each particle in the 
latter case is 5. The fourth and fifth rows contain a loop 
of five particles (moving distances (r^ — Vvif = 5, 1, 2, 1, 
and 1) which extends ('winds') all the way across the lat- 
tice in the x direction. The remaining sites, which share 
the same labels in the left and right, correspond physi- 
cally to particles which have not undergone an exchange. 
The total energy E of the three loops in this configuration 
T'r, is S = (1 + 1) + (5 + 5) (5 + 1 2 + 1 1) = 22. 

The periodic boundary conditions produce an ambi- 
guity in the definition of the energy, since there is more 
than one way in which to compute the distance traveled 
by each particle. We define the energy by choosing the 
smallest of such distances. 

In order to monitor the superfiuid phase transition we 
measure the winding across the lattice in the x and y di- 
rections. Wx counts the difference between the number 
of particles which move to the right across the vertical 
edge of the lattice and those which move to the left. An 
analogous definition applies for Wy across the top, hor- 
izontal edge. We define W"^ = -I- Wy. In the right 
panel of Fig. 1, = 1 and = 0. 

We now introduce our extension of the FK model to 
allow for vacancies. In so doing, we are motivated by 
recent work on 'supersolid' phases^i^ which are, in fact, 
continuations of an extensive history. In the configura- 
tion in the right panel of Fig. 1, each site is labeled by 
one of the 36 sites in the left panel. All sites are occupied. 
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We choose to add the simplest possible term, 



(3lj) (32) (33; (34) (35) (36) 

(25) (26) (27) (23) (29) (30) 

(T9) (20) (2T) (22) (23) (24) 

® ® © © ® ® 

® ® © ® ® 

® ® ® ® ® ® 



(3lj) (32) (33) (34j (35) (36) 

•{30X (26) (27) (23) (22)--(29)— 

(19j (20) 125)— (22)) (23) (24) 

® ® ® ® ® ® 

® ®-® ® ®/® 

® ® ® ® ® ® 



FIG. 1: Left: Unpermuted labeling of the sites in a 6 x 6 lat- 
tice. Right; Representative configuration in our simulation. 
The energy E = 22 and the windings Wx ~ 1, Wy = 0. See 
text for explanation. 
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FIG. 2: Left: Unpermuted labeling of the sites in a 6 x 6 lattice 
when the number of bosons A^;, = 31 < 36. Empty sites are 
labeled by zeroes. Right: Representative configuration in our 
simulation of our extension of the FK model. See text for 
explanation. 



We can define a set of 'vacancies' by assigning '0' to a 
subset of the sites in the lattice (Fig. 2, left). An allowed 
configuration V{Yi) of the system has zeroes at the same 
sites as the vacant sites in r.^ and the occupied site la- 
bels are permuted. An example is given in Fig. 2 (right). 
By a similar calculation to that described for Fig. 1, the 
configuration shown has energy E = 20. Note that we 
will consider "annealed" vacancies: the vacancy density 
is fixed, but the system is allowed to sample all possible 
density locations satisfying that global constraint. 

Once vacancies are allowed, it is possible to introduce 
additional terms in the summand of the partition func- 
tion which control their relative positions on the lattice. 



= +^E(l-^O^rJ(l-5o^r,) . (6) 

(y) 

The sum is over neighboring sites {ij) of the lattice. E.^ 
adds V to the energy for each link which connects sites 
both of which are occupied. As emphasized in Eq. 6, this 
number is the same whether the permuted or unpermuted 
sites are used in the summand, since the collection of 
occupied sites before and after the exchanges is the same. 

The resulting partition function combines the original 
FK exchange term and the new interaction term, 



/?y^(l-,5o.PrJ(l-'5o.Pr, 



(7) 



The inverse temperature (5 appears in the interaction 
term in its usual place. 

We conclude this section by briefiy discussing our sim- 
ulation algorithm. Our approach is a straightforward im- 
plementation of the Metropolis Monte Carlo method^^. 
We suggest a change in our permutation which consists 
of interchanging two, randomly selected, entries Vyi and 
Vvj in the permutation V. If a, vacancy is moved, the list 
of occupied sites and their permutation must be changed 
accordingly. The resulting change A in the argument 
of the exponential appearing in the partition function is 
evaluated, and the change is accepted with probability 
p = min ( 1, ). 

As with most path integral simulations, such local 
moves have difficulty evolving the configuration through 
phase space at large /3 where the important paths are 
dominated by large loops. We therefore also introduce 
'global' moves which shift the elements Vvi by one lat- 
tice constant for sites i across an entire column or row of 
the lattice. Such moves change the vertical or horizontal 
winding of the lattice by AW = ±1. In some simulations 
such moves have low acceptance rates as the system size 
increases. Indeed, this is a primary limitation of simula- 
tions of real Helium^. However, we do not encounter this 
difficulty here. 



III. SIMULATION RESULTS: 
FEYNMAN-KIKUCHI MODEL 

We begin by studying the original FK model. In Fig. 3 
we show data for the mean square winding {W'^) as a 
function of (3 for different lattice sizes. We see that the 
winding becomes nonzero as f3 increases, and that the 
onset of nonzero winding becomes increasingly sharp as 
the lattice size grows. This raw data is suggestive of 
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FIG. 3: Raw data for the mean square winding as a function 
of P on lattice sizes ranging from 8 x 8 to 128 x 128. 
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FIG. 4: Raw data of Fig. 3 for the mean square winding scaled 
by the lattice size. The intersection of the curves determines 
13c = 0.62 ±0.01. 



a critical j3c ~ 0.6 for the development of macroscopic 
loops. In simulations of the FK model in which only local 
exchanges are allowed, on a d = 3 cubic latticeii, Elser 
finds Pc ~ 0.69. Presumably the higher dimensionality 
lowers Pc relative to our d = 2 square lattice while the 
restriction to local exchage would tend to raise f3c- Hence 
a rough match of the critical points is plausible. 

We can make the case for a phase transition, and de- 
termine Pc with higher precision, by scaling our raw data. 
We adopt the usual ansata^i which postulates that the 
dependence of the order parameter on the parameter con- 
trolling the transition and on the lattice size takes the 
scaling form. 



{L,P) = Ly[L\p-P,y 



(8) 



Here / is a universal (lattice size independent) function 
of its argument, and a and b are critical exponents. 
This scaling form is usefully rewritten as. 



(9) 



From this expression it is clear that if we scale the order 
parameter (W^) by the lattice size to an appropriate ex- 
ponent, L~°-{W'^), and plot as a function of the control 
parameter /3, all the curves will cross at the universal 
value /(O) when P — Pc, regardless of the value of the 
second exponent b. 

Fig. 4 presents the results of the analysis in which the 
vertical (order parameter, (W^)) axis alone is scaled. We 
observe a universal crossing of the five curves, and infer 
Pc = 0.62 ± 0.01 and a = 0.61 ± 0.03. 

We calculate the average energy and specific heat from 
the partition function using the standard thermodynam- 
ics formulae. 



d\nZ 
dp 



C 



dE 
df 



(10) 



Fig. [5] shows a plot of specific heat as a function of tem- 
perature. We used the temperature rather than P as our 



horizontal axis, as is more conventional. We see that 
there is a peak in the specific heat at roughly the same 
position as the crossing of the winding Tc = 1/Pc — 
1/0.62 = 1.61. 
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FIG. 5: Plot of specific heat for the original FK model, that 
is, a fully filled lattice with no vacancies and no interactions. 
The value of the critical temperature inferred from the scaling 
of the winding (Fig. 4) is indicated as a vertical dotted line. 
The negative values of C{T) at low T are a finite size artifact, 
as explained in the text. 

The specific heat exhibits a sharp drop, and goes neg- 
ative, close to r = 0. This is a non-physical effect which 
appears due to the finite size of the lattice, as can be seen 
as follows. Starting with the partition fimction of Eq. 5 
we obtain. 



E : 



dlnZ 
dp 



1 EnKe-^ 
P^ Zo 



1 



(E) (11) 



where Zq = X^™ e " is the partition function used in 
Monte Carlo simulations; (E) is our MC calculation of 
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energy. The specific heat 

dE dE 



8(1/(3) 



(12) 



For a finite lattice, as /3 ^ oo, A and B approach con- 
stant values given by all permutations on the Lx L lattice 
having equal probability. As a consequence, as /3 — > cxd, 



1 



1 



1 



/32 /32 



{A-BI3) -^Q- 



So at large /3, or small T, the specific heat becomes neg- 
ative, achieves its minimum, and goes back to zero, as 
T ^ 0. As the lattice size increases, the area of nega- 
tive specific heat moves closer to T = 0. This is because 
the situation when all permutations have approximately 
equal probability occurs when K/f3 ^ 0, or / f3 0, 
so P ^ L^, or T ^ for this area of negative specific 
heat. 

Despite the crudity of the model, we can, following 
Feynman and KikuchiStlfi, use these results to infer a 
rough critical temperature for Helium, The scaling of 
the winding gives f3c = 0.62, or Tc = 1.62. To recover 
physical values for the temperature we note our unit is 
Ti = /ksmcP , where m is the mass of boson, and d 
is lattice spacing. For '^He, using density of 146 kg/rrr', 
and d = 3.57 • 10~i°m, we obtain Ti = 0.95A^ This puts 
our transition temperature around 1.5 K, which is in the 
same ballpark as T\ for Helium. 

We close this section by examining the low T behavior 
of the specific heat, since obtaining the proper exponent 
was the focus of much of the original work on the FK 
model. In three dimensions, where the initial analytic 
studies were performed, a linearly dispersing (phonon) 
mode E{k) ^ ck gives C{T) - at low T. Here 
we are working in d = 2, where instead C(T) ^ T^. 
Fig. [6] shows an attempt to fit C (T) to a power law. The 
least squares fit gives C{T) ~ yi.95±o.i foj. lattices 
of different sizes, which is in reasonable agreement with 
the prediction based on a linearly dispersing mode. (We 
restrict our fit to temperatures higher than those at which 
the finite size artifact negative C(T) values onset.) 



IV. SIMULATION RESULTS: EXTENDED 
FEYNMAN KIKUCHI MODEL 

We now turn to our generalization of the FK model in 
which we allow a lattice with partial filling and a repulsive 
interaction V between nearest neighbor sites. We study 
first the special case of half-filling, where it is possible 
to have perfect ordering of the vacancies/particles in a 
"checkerboard" pattern. 
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FIG. 6: A log-log plot of specific heat versus temperature 
yields data consistent with a straight line of slope 1.95±0.1, in 
agreement with the expected value, 2. We used a temperature 
range 0.5 < T < 1.25, where T is small, but out of range of 
the (non-physical) dip in C, 
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FIG. 7: Raw data for (W^> for the extended FK model with 
half-filling and = 0. 



A. Half-filling: (p) 



How does the introduction of vacancies affect in 
the absence of interactions, I^ = 0? (Henceforth in this 
manuscript we will append a superscipt 'sf to Pc for the 
superfiuid transition, to distinguish it from the P^'^^ for 
charge ordering. See below.) Figs. [7] and [5] are the ana- 
logues of Figs. [3] and m and show the unsealed and scaled 
winding as a function of p. As before, we have done sim- 
ulations for lattices of different sizes to perform the finite 
size scaling. The crossing occurs at Pf = 1.50±0.02. By 
repeating this sweep of P for different V, we can compute 
the superfiuid phase boundary pf{V) in the V — P plane 
at half-filling. This is shown in Fig. [TTl 

A natural question to ask is whether Ey induces va- 
cancy/density ordering. We define the real space density 



6 



20 



15 



10 



8x8 



in 



16x16 - 
32x32 T 
64x64 ♦ 




P 

FIG. 8: Data of Fig. [3 scaled. We infer = 1.50 ± 0.02. 




FIG. 9: The structure factor ^(Tr, tt) is shown as a function 
of /J at half-filling and V = 2.5. 



correlation function, 

c(r,r') = ((p(r)-(p))(p(r')-(p») , (13) 

where the density p(ri) = 1 if the site i is occupied and 
is zero otherwise. The structure function, 

r,r ' 

is the Fourier transform of the density correlations. Here 
N is number of sites. At half-filling, the ordering vector 
q = (7r,7r). In a disordered phase where c(r, r') decays 
exponentially with |r — r'|, S{q) will vanish in the ther- 
modynamic limit (proportional to 1/iV). In an ordered 
phase, S'(q) will go to a constant as the system size in- 
creases, for the appropriate ordering wave vector q. 

Proceeding in analogy with the winding {W'^), we mea- 
sure S{'K,Tr) for different lattice sizes and do finite size 
scaling to determine (3^'^'" for the CDW transition. Rep- 
resentative plots, for V — 2.5, are shown in Figs. [5] and 
[TUl Putting together such sweeps for different interac- 
tion strengths V yields the density order-disorder (CDW) 
phase boundary in the f3 — V plane at half-filling shown 
in Fig. HH 

In fact, a number of aspects of the phase diagram of 
Fig. [11] can be inferred by a mapping to the Ising model. 
We note that the interaction energy Ey of Eq. 6 maps pre- 
cisely to that of the Ising model (the well-known equiva- 
lence of lattice-gas and Ising models) with the definition 
of the spin = 2{p{Ti) — |). The alternating occupied 
and empty site pattern caused by the repulsive particle- 
particle interaction {V > 0) corresponds to an antifer- 
romagnetic arrangement in spin language. The quantity 
V/4: plays the role of the exchange constant J. The criti- 
cal temperature of the Ising model given by the Onsager 
solution Tc — 2.269J then implies a density ordering 
pcdw ^ 4/(2.269 y) = 1.763/y. We expect this resuh 
to be accurate for V » P where Ey dominates over any 
effects on the site densities which might be caused by the 
exchange term. This curve is shown as the dotted line in 
Fig.HH 
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FIG. 10: Scaled data of Fig. [O] The crossing point defines 
Pc'^^ for CDW ordering. The scaling exponent multiplying 
S on the vertical axis is that of the d = 2 Ising model, as 
discussed in the text. 



We can also make a qualitative argument for how the 
phase boundary might bend away from this Ising limit as 
the role of /3 becomes larger. For any given density ar- 
rangement, the exchange of particles provides additional 
configurations of the system associated with permuta- 
tions of the particle indices. Such configurations have 
lower energy when the occupied sites are adjacent. Thus 
we expect that the exchange term will favor "ferromag- 
netic" spin configurations, in competition with the "anti- 
fcrromagnetism" driven by Ey. This suggests a lowering 
of the CDW transition temperature. Just such an in- 
crease in P^*^"" is seen in the phase diagram. One may 
well ask whether there is an extreme limit where the at- 
traction between particles due to the greater ease of ex- 
change becomes so dominant that phase separation oc- 
curs (ferromagnetic clusters in spin language). We will 
address this possibililty later. 

The corresponding effect of V on the superfluid phase 
transition is less easy to describe rigorously, in part be- 
cause we do not begin from a known limit like the On- 
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FIG. 11: The central result of our paper: Phase diagram 
of the extended Feynman-Kikuchi model in the V — (3 plane 
at half-filling. Squares denote simulation results for the su- 
perfluid transition, and triangles the CDW transition. The 
CDW boundary is in good correspondence with the Ising limit 
^cdw ^ 4/(2.269 V), indicated by the dotted line. The circle 
at V = 2.5 is the prediction (3f {p) = Pf(p = l)/p. See text 
for details. 



tions. Since we are interested in the superfluid transi- 
tion, we will restrict ourselves to the case where the two 
particles do exchange, which is reflected in the nonzero 
value of the kinetic energy in the table above. 
The expectation value of the Kinetic energy is 



{K) = 



2-5e 



4e 



2e 



(15) 



Note that (K) = 4 is kinetic energy at V ^ oo. The 
superfluid transition occurs when K ^ (3. If we set (K) = 
f3f and use {K{V ^ oo)) = 4 = (3f {V oo) we find 
that the shift in the superfluid transition is given by, 



PciV ^ oo) 



This result is qualitatively correct aX V ^ 
a small positive shift in (3f relative to V 



(16) 

0, predicting 
oo). 



sager solution in the density transition case. On the one 
hand, increasing V drives the particles apart, making lo- 
cal exchange more expensive, suggesting that (3^^ might 
increase. On the other hand, the superfluid transition 
is not caused by local exchange but instead by global 
winding, and by separating the particles, V might help 
provide regularly spaced "stepping stones," aiding global 
winding and decreasing (3f-. While these qualitative ar- 
guments provide different conclusions, numerically, the 
answer is clear from Fig. [TlJ turning on V decreases (3^^ . 
The effect is not large, however. 

In the limit of large V and half- filling, a perfect CDW 
phase forms, which corresponds to a completely filled 
square lattice with a lattice constant ^/2 larger than the 
original lattice (and rotated by 45 degrees). The squared 
distances in the FK exchange energy will be scaled up by 
a factor of two, and hence (if will be twice the value for 
the no-vacancy FK model of the previous section. Thus 
I3f = 2 (0.62 ± 0.01). This value is shown as the circle 
at = 2.5 in Fig. [Til 

In order to develop a simple understanding of (3^ for 
general V , we consider a small (four site, two particle) 
cluster and enumerate completely the allowed configu- 
rations. There are six possible density configurations, 
which separate into two classes. Four of them have the 
two particles adjacent, while the other two have the two 
particles separated by vacancies: 

Configuration | K energy | P energy | Weight | 
.-.-o-o [1/2 + 9/2 \V |4e"(t+^/^)~ 
• -0-.-0 14/2 + 4/2 |0 I 2e"(^) \ 

We have included the degeneracy factors in the weight. 
For each density configuration, there are two permuta- 



B. Doped system (p) 7^ i 

In this section we consider general filling p. Specifi- 
cally, in Fig. [12] we exhibit the phase diagram in the (3 — p 
plane for two fixed values of the interaction, V — 1.25 and 
V = 2.50. Fig. \n\ was obtained using the same analysis 
as in the earlier sections: Evaluation and scaling of the 
winding and structure factor as a function of (3 for dif- 
ferent V and p. Several features are immediately appar- 
ent from the phase diagram: Charge ordering is, as ex- 
pected, favored close to half-filling, with the highest tran- 
sition temperature at p — 1/2. Interestingly, the shape 
and size of the density ordered region around half-filling 
is in rough agreement with the boundaries obtained for 
checkerboard solid order in the extended boson-Hubbard 
mode l^^i^^ where similar superfluid and charge ordered 
phases are present in an explicitly quantum model. 

In the preceding section we argued that (3^ for half- 
filling and = 00 should be a factor of two larger than 
for the original, no vacancy FK model, and showed this 
was borne out numerically. One might expect that a sim- 
ilar result would be true for general fillings and that (3^ 
would be increased by a factor of 1/p, since this factor 
reflects the increase in the square of the average inter- 
particle spacing. 

/3f(p,y^oo) = (17) 
P 

However, on further consideration, it is not quite so. 
Half-filling and = 00 is a special case and in general 
the particles are not dispersed uniformly, that is, they 
no longer all have the same distance from their nearest 
neighbors. Nevertheless, this relation provides a reason- 
able guide to the density dependence of the superfluid 
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FIG. 12: Phase diagram in the p — jS plane for the extended 
FK model for different values of V. Triangles: the numerically 
obtained cdw transitions. Squares: the numerically obtained 
superfluid transitions. The superfluid phase boundary is rea- 
sonably well approximated by oc 1/p (solid curve). See 
text. 

transition, and is shown on phase diagram at Fig. [T^] as 
the line ''N/SF V = Inf " . In Fig. [H we also see that 
as fi is increased at p = 1/2 and fixed V — 1.25, we go 
from normal to superfluid to supersoUd, where the phases 
are labeled by the behavior of the two order parameters, 
lyV"^) and S{tt,tt). Likewise at larger V — 2.50 we go 
from normal to CDW to supersolid. 

While the order parameters provide unambiguous iden- 
tification of the phases, it is also interesting to see if the 
thermodynamics can pick up the two successive transi- 
tions in the form of separate peaks in the speciflc heat. 
To study C(T), we proceed similarly to the derivation 
for the partition function with only the kinetic energy 
term. Now our partition function has both potential and 
kinetic energy. 

Z = ^e-^'^-f 

n 

dp 




Fig. [m shows the specific heat as a function of tem- 
perature for half-fiUing and interaction V = 1.25. Ac- 
cording to the phase diagram Fig. [T^l the two transition 
points are N/SF at /? = 1.26 (T = 0.79) and N/CDW at 
(3 = 1.47 (T = 0.68). As can be seen, we cannot resolve 
separate peaks in C(T) associated with these transitions. 
It is likely that the critical temperatures are too close and 



.4 
T 

FIG. 13: Specific heat for half-filling, V = 1.25. From Fig.[l2l 
the SF transition is at T = 0.79, and the CDW transition 
at r = 0.68. These values are shown as vertical dashed lines. 
The specific heat cannot resolve the two peaks. 




T 

FIG. 14: Specific heat for half-filling and V = 10. Two peaks 
corresponding to SF and CDW transitions are clearly seen. 
The critical values given by the winding and structure factor 
are shown as vertical dashed lines. The Ising mapping would 
give TcDw = 5.67. 



that the finite size rounding blurs the two peaks into a 
single maximum. 

We can however exhibit separate SF and CDW peaks 
in the specific heat if we push the transitions apart suf- 
ficiently. For example, at V = 10 the CDW transition 
occurs at a much higher temperature than the SF tran- 
sition. Indeed, in Fig. [T3] we can now observe separate 
signatures of the two transitions in C{T). The maxima 
occur close to the transition points given by the order 
parameters. 

Our final results concern the possibility of phase sepa- 
ration. One might argue that in a model with vacancies, 
especially at low or vanishing V , the particles will clump 
together in order to facilitate exchange. Indeed, phase 
separation has been observed in a related model: the 
bose Hubbard Hamiltonian with ring exchange, precisely 
due to this mechanism^. 
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FIG. 15: ^(gO for various V, L = 16, (3 = 1.5. S'(27r/L,0) 
should signal phase separation if it is ~ 0.1. It is however 
not close to that number for V > 0. For comparison, plots 
for V < are given, and they indicate phase separation. 



Phase separation is signalled by a peak in the density 
structure factor at small momenta q (as opposed to the 
CDW ordering vector at the largest q = (tt, tt)). Crudely 
speaking there are real space density fluctuations at long 
wavelengths, corresponding to a lattice with one side half 
occupied and the other half empty. These translate into 
a peak in S{q) at small q. Note that in a canonical en- 
semble simulation such as is performed here we cannot 
set q = (0, 0) since that value of the structure factor is 
just a constant set by the filling. Indeed with our defi- 
nition of the density correlations in terms of fluctuations 
about the average density per site, Eq. 14, 5'(0,0) = 0. 

For a perfectly phase separated state with all particles 
on the right-most half of the lattice we find S{2tt/L, 0) w 
0.11 for L = 16. In Fig. [15] we see that S{2Tr/L,0) is 
^ 1/100 of that figure for > 0. We conclude there is 
no phase separation in this model. The absence of a sig- 
nal for phase separation is in contrast to the behavior in 
the related bose-Hubbard model with ring exchange^^-^* 
There the same quantity, the average of the structure fac- 



tors at the three lowest momenta shows a sharp rise with 
increased exchange, as one enters the superfluid phase. 



Conclusions 

In this paper we have presented Monte Carlo simula- 
tions of the phase diagram of an extension of the d = 2 
Feynman-Kikuchi model which includes vacancies. We 
found phases which have density and superfluid order, 
and where these two types of order coexist. Unlike the 
boson-Hubbard model where supersolid order requires 
doping away from half-filling in the extended FK model 
Ps is nonzero even in the defect-free checkerboard solid. 
The reason is that the boson-Hubbard kinetic energy 
moves particles only between near-neighbor sites. Bosons 
cannot exchange without passing through an energeti- 
cally unfavorable region. But in the FK model, exchange 
at longer range can occur without ever 'passing through' 
the rare configurations with near-neighbor sites that are 
occupied. One might expect that in the FK model which 
is restricted to local exchange, the half-filled supersolid 
might be eliminated. 

A further problem of interest in the extended FK is 
to consider "quenched" vacancies in which the locations 
of the empty sites are frozen throughout the simulation. 
Here again we might expect that when a restriction to 
local exchange is enforced, there could be a destruction 
of the superfluid transition as the percolation threshold 
is crossed. In the model allowing exchanges of arbitrary 
distance, one expects a more trivial increase in /3c, but 
that the superfluid transition would likely persist. 

A final avenue for exploration would be the inclusion 
of a one-body vacancy potential which could be chosen to 
confine the particles preferentially towards the center of 
the lattice. Such simulations would connect with recent 
experiments on cold atoms in magnetic and laser traps. 

We acknowledge support from the National Science 
Foundation under award NSF ITR 0313390, and useful 
input from T.Tremeloes. 
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